from __future__ import division # use so 1/2 = 0.5, etc.
%pylab inline
import sk_dsp_comm.sigsys as ss
import scipy.signal as signal
from IPython.display import Image, SVG

Populating the interactive namespace from numpy and matplotlib

pylab.rcParams['savefig.dpi'] = 100 # default 72
%config InlineBackend.figure_formats=['svg'] # SVG inline viewing

Introduction to Python and the Jupyter Notebook

t = arange(-4,4,.01)
x = cos(2*pi*t)
plot(t,x)

grid()

Rectangle and Triangle Pulses Defined

Before showing more examples, consider some familiar signal primitives in your signals and systems background.

t = arange(-5,5,.01)
x_rect = ss.rect(t-3,2)
x_tri = ss.tri(t+2,1.5)
subplot(211)
plot(t,x_rect)
grid()
ylabel(r'$\Pi((t-3)/2)$');
subplot(212)
plot(t,x_tri)
grid()
xlabel(r'Time (s)')
ylabel(r'$\Lambda((t+2)/1.5)$');
tight_layout()

• Consider an interactive version of the above:

# Make an interactive version of the above
from ipywidgets import interact, interactive

def pulses_plot(D1,D2,W1,W2):
    t = arange(-5,5,.01)
    x_rect = ss.rect(t-D1,W1)
    x_tri = ss.tri(t-D2,W2)
    subplot(211)
    plot(t,x_rect)
    grid()
    ylabel(r'$\Pi((t-3)/2)$');
    subplot(212)
    plot(t,x_tri)
    grid()
    xlabel(r'Time (s)')
    ylabel(r'$\Lambda((t+2)/1.5)$');
    tight_layout()

interactive_plot = interactive(pulses_plot,D1 = (-3,3,.5), D2 = (-3,3,.5), W1 = (0.5,2,.25), W2 = (0.5,2,.25));
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot

More Signal Plotting

The basic pulse shapes (primitives) defined in the module ssd.py are very useful for working Text 2.13a &d, but there are also times when you need a custom piecewise function.

Simple Cases:

Consider plotting

• \(x_1(t) = \sin(2\pi\cdot 5t) \Pi((t-2)/2)\) for \(0\leq t \leq 10\)
• \(x_2(t) = \sum_{n=-\infty}^\infty = \Pi((t-5n)/1)\) for \(-10 \leq t \leq 10\)

t1 = arange(0,10+.01,.01) # arange stops one step size less than the upper limit
x1 = sin(2*pi*5*t1)* ss.rect(t1-2,2)
subplot(211)
plot(t1,x1)
xlabel(r'Time (s)')
ylabel(r'$x_1(t)$')
grid()
t2 = arange(-10,10,.01)
# Tweak mod() to take on negative values
x2 = ss.rect(mod(t2+2.5,5)-2.5,1)
subplot(212)
plot(t2,x2)
xlabel(r'Time (s)')
ylabel(r'$x_2(t)$')
grid()
tight_layout()

def x3_func(t):
    """
    Create a piecewise function for plotting x3
    """
    x3 = zeros_like(t)
    for k,tk in enumerate(t):
        if tk >= 0 and tk <= 3:
            x3[k] = 1 + tk**2
        elif tk > 3 and tk <= 5:
            x3[k] = cos(2*pi*5*tk)
    return x3

t3 = arange(-2,6+.01,.01)
x3 = x3_func(t3)
plot(t3,x3)
xlabel(r'Time (s)')
ylabel(r'$x_3(t)$')
xlim([-2,6])
grid()

26/2

13.0

Energy and Power Signals

You can create an approximation to the waveform over a finite number of periods by doing a little programming:

def periodic_tri(t,tau,T,N):
    """
    Approximate x2(t) by running the sum index from -N to +N.
    The period is set by T and tau is the tri pulse width
    parameter (base width is 2*tau).

    Mark Wickert January 2015
    """
    x = zeros_like(t)
    for n in arange(-N,N+1):
        x += ss.tri(t-T*n,tau)
    return x

t = arange(-10,10,.001)
x = periodic_tri(t,2,6,10)
plot(t,x)
plot(t,abs(x)**2)
grid()
#xlim([-5,5])
xlabel(r'Time (s)')
ylabel(r'$x_2(t)$ and $x_2^2(t)$');

For the power calculation create a time array that runs over exactly one period. Below is the case for the original problem statement.

T0 = 6
tp = arange(-T0/2,T0/2+.001,.001)
xp = periodic_tri(tp,2,T0,5)
plot(tp,xp)
plot(tp,abs(xp)**2)
legend((r'$x(t)$', r'$|x(t)|^2$'),loc='best',shadow=True)
grid();
xlim([-T0/2,T0/2])
xlabel(r'Time (s)')
ylabel(r'$x_2(t)$ and $x_2^2(t)$');

#Power calculation
Px2 = (1/T0)*sum(xp**2)*.001 # rectangular partitions for integral
print('Power estimate via numerical integration: %2.4f W' % Px2)

Power estimate via numerical integration: 0.2222 W

Power in the Sum of Two Sinusoids

t = arange(-10,10,.001)
x1 = 4*cos(2*pi*10*t)
x2 = 3*cos(2*pi*3.45*t+pi/9)
plot(t,x1)
plot(t,x2)
plot(t,x1+x2)
grid()
xlabel(r'Time (s)')
ylabel(r'Amplitude')
legend((r'$x_1(t)$', r'$x_2(t)$', r'$x_1(t)+x_2(t)$'),loc='best',shadow=True)
xlim([-.1,.1]);

print('Power calculations: %3.2f, %3.2f, %3.2f' \
      % (var(x1),var(x2),var(x1+x2)))

Power calculations: 8.00, 4.50, 12.50

print('Theory: %3.2f, %3.2f, %3.2f' \
      % (4**2/2,3**2/2,4**2/2+3**2/2))

Theory: 8.00, 4.50, 12.50

Fourier Series and Line Spectra Plotting

Being able to easily plot the line spectra of periodic signals will hopefully enhance your understanding. The module ss.py contains the function ss.line_spectra() for this purpose. The function assumes that the Fourier coefficients, \(X_n\) are available for a real signal \(x(t)\). The function plots line spectra as: * The two-sided magnitude spectra * The two-sided magnitude spectra in dB with an adjustable floor level in dB * The two-sided phase spectra in radians * The one-sided line spectra corresponding to the three cases listed immediately above Examples are given below for the case of a simple pulse train and then for a trapezoidal pulse train. IN the case of the trapezoidal pulse train the underlying Fourier coefficients are obtained numerically using the FFT as described in the course notes.

A fundamental requirement in using ss.line_spectra() is to beable to supply the coefficients starting with the DC term coefficient \(X_0\) and moving up to the \(N\)th harmonic. Before plotting the pulse train line spectra I first describe a helper function for visualizing the pulse train waveform.

Pulse Train

def pulse_train(Np,fs,tau,t0):
    """
    Generate a discrete-time approximation to a continuous-time
    pulse train signal. Amplitude values are [0,1]. Scale and offset
    later if needed.

    Inputs
    ------
     Np = number of periods to generate
     fs = samples per period
    tau = duty cycle
     t0 = pulse delay time relative to first rising edge at t = 0

    Return
    ------
    t = time axis array
    x = waveform

    Mark Wickert, January 2015
    """
    t = arange(0,Np*fs+1,1)/fs #time is normalized to make period T0 = 1.0
    x = zeros_like(t)
    # Using a brute force approach, just fill x with the sample values
    for k,tk in enumerate(t):
        if mod(tk-t0,1) <= tau and mod(tk-t0,1) >= 0:
            x[k] = 1
    return t,x

tau = 1/8; fs = 8*16; t0 = 0 # note t0 = tau/2
subplot(211)
t,x = pulse_train(4,fs,tau,t0)
plot(t,x) # Just a plot of xa(t)
ylim([-.1,1.1])
grid()
ylabel(r'$x_a(t)$')
title(r'Pulse Train Signal: (top) $x_a(t)$, (bot) $x_b(t) = 1-x_a(t)$');
subplot(212)
t,x = pulse_train(4,fs,tau,t0)
plot(t,1-x) # Note here y(t) = 1 - x(t), a special case of
ylim([-.1,1.1]) # y(t) = A + B*x(t) in the notes
grid()
xlabel(r'Time ($t/T_0$)')
ylabel(r'$x_b(t)$');

Example: Pulse Train Line Spectra

n = arange(0,25+1) # Get 0 through 25 harmonics
tau = 0.125; f0 = 1; A = 1;
Xn = A*tau*f0*sinc(n*f0*tau)*exp(-1j*2*pi*n*f0*tau/2)
# Xn = -Xn # Convert the coefficients from xa(t) t0 xb(t)
# Xn[0] += 1
figure(figsize=(6,2))
f = n # Assume a fundamental frequency of 1 Hz so f = n
ss.line_spectra(f,Xn,mode='mag',sides=2,fsize=(6,2))
xlim([-25,25]);
#ylim([-50,10])
figure(figsize=(6,2))
ss.line_spectra(f,Xn,mode='phase',fsize=(6,2))
xlim([-25,25]);

<Figure size 432x144 with 0 Axes>

<Figure size 432x144 with 0 Axes>

Example: Trapezoidal Pulse

Consider the line spectra of a finite rise and fall time pulse train is of practical interest. The function trap_pulse() allows you first visualize one period of the trapezoidal pulse train, and then use this waveform in obtaining numerically the Fourier coefficients of this signal. PLotting the corresponding line spectra follows.

A point to be main is that by slowing down the edges (rise time/fall time) of the pulse train the amplitude of the harmonics falls off more rapidly. When considering the clock speed in todays PCs this can be a good thing as harmonic emission is an issue.

def trap_pulse(N,tau,tr):
    """
    xp = trap_pulse(N,tau,tr)

    Mark Wickert, January 2015
    """
    n = arange(0,N)
    t = n/N
    xp = zeros(len(t))
    # Assume tr and tf are equal
    T1 = tau + tr
    # Create one period of the trapezoidal pulse waveform
    for k in n:
        if t[k] <= tr:
            xp[k] = t[k]/tr
        elif (t[k] > tr and t[k] <= tau):
            xp[k] = 1
        elif (t[k] > tau and t[k] < T1):
            xp[k] = -t[k]/tr + 1 + tau/tr;
        else:
            xp[k] = 0
    return xp, t

Let \(\tau = 1/8\) and \(t_r = 1/20\):

# tau = 1/8, tr = 1/20
N = 1024
xp,t = trap_pulse(N,1/8,1/20)
Xp = fft.fft(xp)
figure(figsize=(6,2))
plot(t,xp)
grid()
title(r'Spectra of Finite Risetime Pulse Train: $\tau = 1/8$ $t_r = 1/20$')
ylabel(r'$x(t)$')
xlabel('Time (s)')
f = arange(0,N/2)
ss.line_spectra(f[0:25],Xp[0:25]/N,'magdB',floor_dB=-80,fsize=(6,2))
ylabel(r'$|X_n| = |X(f_n)|$ (dB)');
#% tau = 1/8, tr = 1/10
xp,t = trap_pulse(N,1/8,1/10)
Xp = fft.fft(xp)
figure(figsize=(6,2))
plot(t,xp)
grid()
title(r'Spectra of Finite Risetime Pulse Train: $\tau = 1/8$ $t_r = 1/10$')
ylabel(r'$x(t)$')
xlabel('Time (s)')
ss.line_spectra(f[0:25],Xp[0:25]/N,'magdB',floor_dB=-80,fsize=(6,2))
ylabel(r'$|X_n| = |X(f_n)|$ (dB)');

With the edge speed slowed down it is clear that the harmonics drop off faster.

Fourier Transforms

def FT_approx(x,t,Nfft):
    '''
    Approximate the Fourier transform of a finite duration
    signal using scipy.signal.freqz()

    Inputs
    ------
       x = input signal array
       t = time array used to create x(t)
    Nfft = the number of frdquency domain points used to
           approximate X(f) on the interval [fs/2,fs/2], where
           fs = 1/Dt. Dt being the time spacing in array t

    Return
    ------
    f = frequency axis array in Hz
    X = the Fourier transform approximation (complex)

    Mark Wickert, January 2015
    '''
    fs = 1/(t[1] - t[0])
    t0 = (t[-1]+t[0])/2 # time delay at center
    N0 = len(t)/2 # FFT center in samples
    f = arange(-1/2,1/2,1/Nfft)
    w, X = signal.freqz(x,1,2*pi*f)
    X /= fs # account for dt = 1/fs in integral
    X *= exp(-1j*2*pi*f*fs*t0)# time interval correction
    X *= exp(1j*2*pi*f*N0)# FFT time interval is [0,Nfft-1]
    F = f*fs
    return F, X

Example: Rectangular Pulse

fs = 100 # sampling rate in Hz
tau = 1
t = arange(-5,5,1/fs)
x0 = ss.rect(t-.5,tau)
figure(figsize=(6,5))
subplot(311)
plot(t,x0)
grid()
ylim([-0.1,1.1])
xlim([-2,2])
title(r'Exact Waveform')
xlabel(r'Time (s)')
ylabel(r'$x_0(t)$');

# FT Exact Plot
fe = arange(-10,10,.01)
X0e = tau*sinc(fe*tau)
subplot(312)
plot(fe,abs(X0e))
#plot(f,angle(X0))
grid()
xlim([-10,10])
title(r'Exact Spectrum Magnitude')
xlabel(r'Frequency (Hz)')
ylabel(r'$|X_0e(f)|$');

# FT Approximation Plot
f,X0 = ss.ft_approx(x0,t,4096)
subplot(313)
plot(f,abs(X0))
#plot(f,angle(X0))
grid()
xlim([-10,10])
title(r'Approximation Spectrum Magnitude')
xlabel(r'Frequency (Hz)')
ylabel(r'$|X_0(f)|$');
tight_layout()

Example: Text Problem 2.31a Drill Down

fs = 100
t = arange(-5,5,1/fs)
x1 = ss.rect(t+1/2,1)-ss.rect(t-1/2,1)
subplot(211)
plot(t,x1)
grid()
ylim([-1.1,1.1])
xlim([-2,2])
xlabel(r'Time (s)')
ylabel(r'$x_1(t)$');
fe = arange(-10,10,.01)
X1e = 2*1j*sinc(fe)*sin(pi*fe)
f,X1 = ss.ft_approx(x1,t,4096)
subplot(212)
plot(f,abs(X1))
plot(fe,abs(X1e))
#plot(f,angle(X1))
legend((r'Num Approx',r'Exact'),loc='best')
grid()
xlim([-10,10])
xlabel(r'Frequency (Hz)')
ylabel(r'$|X_1(f)|$');
tight_layout()

• Notice the numerical approximation and exact spectral plots overlay one another

Example: Modulation Theorem

fs = 100 # sampling rate in Hz
tau = 1
t = arange(-5,5,1/fs)
x3 = ss.tri(t,tau)
y = x3*cos(2*pi*10*t)
subplot(211)
plot(t,x3)
plot(t,y)
grid()
ylim([-1.1,1.1])
xlim([-2,2])
legend((r'$x_3(t)$', r'$y(t)$'),loc='lower right',shadow=True)
title(r'Time Domain: $x_3(t)$ and $y(t)=x_3(t)\cos(2\pi\cdot 5\cdot t)$')
xlabel(r'Time (s)')
ylabel(r'$y(t)$');
f,Y = ss.ft_approx(y,t,4096)
subplot(212)
plot(f,abs(Y))
#plot(f,angle(X0))
grid()
title(r'Frequency Domain: $Y(f)$')
xlim([-15,15])
xlabel(r'Frequency (Hz)')
ylabel(r'$|Y(f)|$');
tight_layout()

Example: Representing a Bandlimited Signal

fs = 100 # sampling rate in Hz
W = 5
t = arange(-5,5,1/fs)
x4 = 2*W*sinc(2*W*t)
figure(figsize=(6,2))
plot(t,x4)
grid()
#ylim([-1.1,1.1])
xlim([-2,2])
title(r'Time Domain: $x_4(t),\ W = 5$ Hz')
xlabel(r'Time (s)')
ylabel(r'$x_4(t)$');
f,X4 = ss.ft_approx(x4,t,4096)
figure(figsize=(6,2))
plot(f,abs(X4))
grid()
title(r'Frequency Domain: $X_4(f)$')
xlim([-10,10])
xlabel(r'Frequency (Hz)')
ylabel(r'$|X_4(f)|$');
figure(figsize=(6,2))
plot(f,20*log10(abs(X4)))
grid()
title(r'Frequency Domain: $X_4(f)$ in dB')
ylim([-50,5])
xlim([-10,10])
xlabel(r'Frequency (Hz)')
ylabel(r'$|X_4(f)|$ (dB)');

Note: The dB version (last plot) reveals that the first sidelobes of the spectrum are only down ~21dB. Increasing the length of the time window will not help. The spectral side lobes will become more tightly packed, but the first sidelobe will still be down only 21dB. With other pulse shapes in the time domain, i.e., note simply a truncted \(\text{sinc}()\) function reduced sidelobes can be obtained.

Convolution

t = arange(-2,2.001,.001)
p1 = ss.rect(t,1)
p2 = ss.rect(t,3)
y,ty = ss.conv_integral(p1,t,p2,t)
plot(ty,y)
ylim([-.01,1.01])
grid()
xlabel(r'Time (s)')
ylabel(r'$y(t)$');

Output support: (-4.00, +4.00)

For convolutions involving semi-infinite signals, such as \(u(t)\), you can tell ssd.conv_integral() about this via the optional extent argument. See the function help using

ss.conv_integral?

# Consider a pulse convolved with an exponential ('r' type extent)
tx = arange(-1,8,.01)
x = ss.rect(tx-2,4) # pulse starts at t = 0
h = 4*exp(-4*tx)*ss.step(tx)
y,ty = ss.conv_integral(x,tx,h,tx,extent=('f','r')) # note extents set
plot(ty,y) # expect a pulse charge and discharge waveform
grid()
title(r'$\Pi((t-2)/4)\ast 4 e^{-4t} u(t)$')
xlabel(r'Time (s)')
ylabel(r'$y(t)$');

Output support: (-2.00, +6.99)

Spectrum of PN Sequence (exact)

The cell below is a copy of the earlier pulse train line spectra example. Use this as a template to create the solution to the PN code problem of HW 3.

n = arange(0,25+1) # Get 0 through 25 harmonics
tau = 0.125; f0 = 1; A = 1;
Xn = A*tau*f0*sinc(n*f0*tau)*exp(-1j*2*pi*n*f0*tau/2)
# Xn = -Xn # Convert the coefficients from xa(t) t0 xb(t)
# Xn[0] += 1
figure(figsize=(6,2))
f = n # Assume a fundamental frequency of 1 Hz so f = n
ss.line_spectra(f,Xn,mode='mag',sides=2,fsize=(6,2))
xlim([-25,25]);
#ylim([-50,10])
figure(figsize=(6,2))
ss.line_spectra(f,Xn,mode='phase',fsize=(6,2))
xlim([-25,25]);

<Figure size 432x144 with 0 Axes>

<Figure size 432x144 with 0 Axes>

Spectrum of PN Sequence (approx)

The code below approximates the PSD of the PN code using a numerical approximation to the Fourier coefficients, \(X_n\). This development may be useful for the lab, as you can esily change the waveform level without having to rework the theory.

The approach taken here to create one period of the PN waveform at 10 samples per bit. The line containing the function ss.upsample() converts the bit sequence into a waveform by upsampling and filtering with a rectangular pulse shape (ones(10)). The function ss.fs_coeff() numerically calculates the \(X_n\)‘s. To plot the PSD from the Fourier coefficients we use

\[ \begin{align}\begin{aligned} S_x(f) = \sum_{n=-\infty}^\infty |X_n|^2 \delta(f-nf_0)\\where :math:`f_0` in this case is :math:`1/(MT_0)` with :math:`T_0` beging the bit period and :math:`M` the code period in bits.\end{aligned}\end{align} \]

x_PN4 = ss.m_seq(4)
x = signal.lfilter(ones(10),1,ss.upsample(x_PN4,10))
t = arange(0,len(x))/10
figure(figsize=(6,2))
plot(t,x);
title(r'Time Domain and PSD of $M=15$ PN Code with $T = 1$')
xlabel(r'Time (s)')
ylabel(r'x(t)')
axis([0,15,-0.1,1.1]);
grid()
# 10 samples/bit so 150 samples/period
# harmonics spaced by 1/(15*T) = 1/15
Xk,fk = ss.fs_coeff(x,45,1/15)
ss.line_spectra(fk,Xk,'magdB',lwidth=2.0,floor_dB=-50,fsize=(6,2))
xlim([-3,3])
ylabel(r'$|X_n| = |X(f_n)|$ (dB)');

# Line spacing
1/15

0.06666666666666667

import sk_dsp_comm.digitalcom as dc
y_PN5_bits = ss.PN_gen(10000,5)
# Convert to waveform level shifted to +/-1 amplitude
y = 2*signal.lfilter(ones(10),1,ss.upsample(y_PN5_bits,10))-1
# Find the time averaged autocorrelation function normalized
# to have a peak amplitude of 1
Ry,tau = dc.xcorr(y,y,400)
# We know Ry is real so strip small imag parts from FFT-based calc
Ry = Ry.real

fs = 10
t = arange(len(y))/fs
plot(t[:500],y[:500])
title(r'PN Waveform for 5 Stages (Period $2^5 -1 = 31$ bits)')
ylabel(r'Amplitude')
xlabel(r'Bits (10 samples/bit)')
grid();

tau_s = tau/10
figure(figsize=(6,2))
plot(tau_s,Ry)
title(r'Autocorrelation and PSD Estimates for $M=31$ with $T = 1$')
xlabel(r'Autocorrelation Lag $\tau$ (s)')
ylabel(r'$R_y(\tau)$')
grid();
figure(figsize=(6,2))
psd(y,2**12,10)
xlabel(r'Frequency (Hz)')
ylabel(r'$S_y(f)$ (dB)')
#xlim([0,.002]);
ylim([-30,20]);

In Lab 2 of ECE 4670 a C/C++ version of a PN generator is implemented to run the ARM mbed LPC 1768 microcontroller (https://www.sparkfun.com/products/9564). At the heart of this code is:

// Globals defined as unsigned int
tap1 -= 1;
tap2 -= 1;
mask1 = 0x1 << (tap1);
mask2 = 0x1 << (tap2);
bit = 0x0;
sync = 0x0;

void gen_PN() {
    my_pin5 = bit;
    my_pin6 = synch_bit;
    led2 = bit;
    led3 = synch_bit;
    if (clk_state == 0x1)
    {
    // Advance m-sequence generator by one bit
    // XOR tap1 and tap2 SR values and feedback to input
        fb = ((sr & mask1)>> tap1) ^ ((sr & mask2) >> tap2);
        sr = (sr << 1) + fb;
        bit = sr & 0x1;
        // Use random number generator in place of m-sequence bits
        if (DIP_sw4)
        {
            bit = rand_int() & 0x1;
        }
        clk_state = 0x0;
        // See if all 1's condition exists in SR
        if ((sr & synch) == synch) {
            synch_bit = 0x1;
        }
        else
        {
            synch_bit = 0x0;
        }
    }
    else
    {
        if (DIP_sw1) bit = !bit;
        clk_state = 0x1;
    }
}

The data type is unsigned int, which on the mbed is uint16_t, an unsigned 16-bit integer. A single unsigned integer is used as a 16-bit shift register with the LSB, furthest bit to the right, used to represent the first register stage. The shift register is advanced using a left shift << bitwise operation. We can code this Python almost directly, as shown below.

class bitwise_PN(object):
    """
    Implement a PN generator using bitwise manipulation for
    the shift register. The LSB holds b0 and bits are shifted left.
              +----+----+----+----+----+----+----+
         sr = |bM-1| .. |bM-k| .. | b2 | b1 | b0 |
              +----+----+----+----+----+----+----+
                 |         |
    Feedback:(tap1-1)   (tap2-1)    Shift left using <<

    Mark Wickert February 2017
    """
    def __init__(self,tap1,tap2,Nstage,sr_initialize):
        """
        Initialize the PN generator object
        """
        self.tap1 = tap1 - 1
        self.tap2 = tap2 - 1
        self.mask1 = 0x1 << (tap1 - 1) # to select bit of interest
        self.mask2 = 0x1 << (tap2 - 1) # to select bit of interest
        self.Nstage = Nstage
        self.period = 2**Nstage - 1
        self.sr = sr_initialize
        self.bit = 0
        self.sync_bit = 0

    def clock_PN(self):
        '''
        Method to advance m-sequence generator by one bit
        XOR tap1 and tap2 SR values and feedback to input
        '''
        fb = ((self.sr & self.mask1)>> self.tap1) ^ \
            ((self.sr & self.mask2) >> self.tap2)
        self.sr = (self.sr << 1) + fb
        self.sr = self.sr & self.period # set MSBs > Nstage to 0
        self.bit = self.sr & 0x1 # output LSB from SR
        # See if all 1's condition exits in SR, if so output a synch pulse
        if ((self.sr & self.period) == self.period):
            self.sync_bit = 0x1
        else:
            self.sync_bit = 0x0
        print('output = %d, sr contents = %s, sync bit = %d' \
             % (self.bit, binary(self.sr, self.Nstage), self.sync_bit))

# A simple binary format display function which shows
# leading zeros to a fixed bit width
def binary(num, length=8):
    return format(num, '#0{}b'.format(length + 2))

PN1 = bitwise_PN(10,7,10,0x1)

PN1.clock_PN()

output = 0, sr contents = 0b0000000010, sync bit = 0

# sr initial condition
sr = 0b1

Nout = 20
x_out = zeros(Nout)
s_out = zeros(Nout)
PN1 = bitwise_PN(3,2,3,0x1)
for k in range(Nout):
    PN1.clock_PN()
    x_out[k] = PN1.bit
    s_out[k] = PN1.sync_bit

output = 0, sr contents = 0b010, sync bit = 0
output = 1, sr contents = 0b101, sync bit = 0
output = 1, sr contents = 0b011, sync bit = 0
output = 1, sr contents = 0b111, sync bit = 1
output = 0, sr contents = 0b110, sync bit = 0
output = 0, sr contents = 0b100, sync bit = 0
output = 1, sr contents = 0b001, sync bit = 0
output = 0, sr contents = 0b010, sync bit = 0
output = 1, sr contents = 0b101, sync bit = 0
output = 1, sr contents = 0b011, sync bit = 0
output = 1, sr contents = 0b111, sync bit = 1
output = 0, sr contents = 0b110, sync bit = 0
output = 0, sr contents = 0b100, sync bit = 0
output = 1, sr contents = 0b001, sync bit = 0
output = 0, sr contents = 0b010, sync bit = 0
output = 1, sr contents = 0b101, sync bit = 0
output = 1, sr contents = 0b011, sync bit = 0
output = 1, sr contents = 0b111, sync bit = 1
output = 0, sr contents = 0b110, sync bit = 0
output = 0, sr contents = 0b100, sync bit = 0

stem(x_out)
stem(0.2*s_out,markerfmt = 'ro')
ylim([0,1.1])

(0, 1.1)

Cross Correlation and Signal Delay

import sk_dsp_comm.digitalcom as dc
x_PN4_bits = ss.PN_gen(10000,6)
# Convert to waveform level shifted to +/-1 amplitude
x_s = 2*signal.lfilter(ones(10),1,ss.upsample(x_PN4_bits,10))-1
# Form a delayed version of x_S
T_D = 35 # 35 sample delay
y_s = signal.lfilter(concatenate((zeros(T_D),array([1]))),1,x_s)
figure(figsize=(6,2))
plot(x_s[:200])
plot(y_s[:200])
ylim([-1.1,1.1])
title(r'Delayed and Undelayed Signals for $T_D = 35$ Samples')
xlabel(r'Samples (10/PN bit)')
ylabel(r'$x_s(t)$ and $y_s(t)$')
grid();
# Find the time averaged autocorrelation function normalized
# to have a peak amplitude of 1
Ryx,tau = dc.xcorr(y_s,x_s,200) #note order change
# We know Ryx is real
Ryx = Ryx.real
tau_s = tau/10
figure(figsize=(6,2))
plot(tau_s,Ryx)
title(r'Cross Correlation for $M=4$ with $T = 1$ and Delay 35 Samples')
xlabel(r'Autocorrelation Lag $\tau$ (s)')
ylabel(r'$R_{yx}(\tau)$')
grid();

Spectral Containment Bandwidth (text problem 2.55)

Example:

fn = arange(0,10,.001)
Gn = 4*sinc(fn)**2 * sin(pi*fn)**2
Gn_cumsum = cumsum(Gn)
Gn_tot = sum(Gn)
plot(fn,Gn_cumsum/Gn_tot)
grid()
xlabel('Normalized Frequency $f\tau$')
ylabel('Fractional Power Containment');

fn_idx = find(abs(Gn_cumsum/Gn_tot - 0.9)< 0.0005)
fn_idx

/home/docs/checkouts/readthedocs.org/user_builds/scikit-dsp-comm/envs/develop/lib/python3.5/site-packages/ipykernel_launcher.py:1: MatplotlibDeprecationWarning: The find function was deprecated in Matplotlib 2.2 and will be removed in 3.1.
  """Entry point for launching an IPython kernel.

array([1446, 1447, 1448, 1449, 1450])

print('The normalized 90 percent containment bandwidth is %2.2f Hz-s.' \
      % fn[1448])

The normalized 90 percent containment bandwidth is 1.45 Hz-s.

Filter Analysis

To facilitate the performance analysis of both discrete-time and continuous-time filters, the functions freqz_resp() and freqs_resp() are available (definitions below, respectively). With these functions you can quickly move from z-domain or s-domain rational system function coefficients to visualization of the filter frequency response * Magnitude * Magnitude in dB * Phase in radians * Group delay in samples or seconds (digital filter) * Group delay in seconds (analog filter)

def freqz_resp(b,a=[1],mode = 'dB',fs=1.0,Npts = 1024,fsize=(6,4)):
    """
    A method for displaying digital filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freq_resp(self,mode = 'dB',Npts = 1024)

    A method for displaying the filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freqs_resp(b,a=[1],Dmin=1,Dmax=5,mode = 'dB',Npts = 1024,fsize=(6,4))

        b = ndarray of numerator coefficients
        a = ndarray of denominator coefficents
     Dmin = start frequency as 10**Dmin
     Dmax = stop frequency as 10**Dmax
     mode = display mode: 'dB' magnitude, 'phase' in radians, or
            'groupdelay_s' in samples and 'groupdelay_t' in sec,
            all versus frequency in Hz
     Npts = number of points to plot; defult is 1024
    fsize = figure size; defult is (6,4) inches

    Mark Wickert, January 2015
    """
    f = np.arange(0,Npts)/(2.0*Npts)
    w,H = signal.freqz(b,a,2*np.pi*f)
    plt.figure(figsize=fsize)
    if mode.lower() == 'db':
        plt.plot(f*fs,20*np.log10(np.abs(H)))
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Gain (dB)')
        plt.title('Frequency Response - Magnitude')

    elif mode.lower() == 'phase':
        plt.plot(f*fs,np.angle(H))
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Phase (rad)')
        plt.title('Frequency Response - Phase')

    elif (mode.lower() == 'groupdelay_s') or (mode.lower() == 'groupdelay_t'):
        """
        Notes
        -----

        Since this calculation involves finding the derivative of the
        phase response, care must be taken at phase wrapping points
        and when the phase jumps by +/-pi, which occurs when the
        amplitude response changes sign. Since the amplitude response
        is zero when the sign changes, the jumps do not alter the group
        delay results.
        """
        theta = np.unwrap(np.angle(H))
        # Since theta for an FIR filter is likely to have many pi phase
        # jumps too, we unwrap a second time 2*theta and divide by 2
        theta2 = np.unwrap(2*theta)/2.
        theta_dif = np.diff(theta2)
        f_diff = np.diff(f)
        Tg = -np.diff(theta2)/np.diff(w)
        max_Tg = np.max(Tg)
        #print(max_Tg)
        if mode.lower() == 'groupdelay_t':
            max_Tg /= fs
            plt.plot(f[:-1]*fs,Tg/fs)
            plt.ylim([0,1.2*max_Tg])
        else:
            plt.plot(f[:-1]*fs,Tg)
            plt.ylim([0,1.2*max_Tg])
        plt.xlabel('Frequency (Hz)')
        if mode.lower() == 'groupdelay_t':
            plt.ylabel('Group Delay (s)')
        else:
            plt.ylabel('Group Delay (samples)')
        plt.title('Frequency Response - Group Delay')
    else:
        s1 = 'Error, mode must be "dB", "phase, '
        s2 = '"groupdelay_s", or "groupdelay_t"'
        print(s1 + s2)

def freqs_resp(b,a=[1],Dmin=1,Dmax=5,mode = 'dB',Npts = 1024,fsize=(6,4)):
    """
    A method for displaying analog filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freqs_resp(b,a=[1],Dmin=1,Dmax=5,mode='dB',Npts=1024,fsize=(6,4))

        b = ndarray of numerator coefficients
        a = ndarray of denominator coefficents
     Dmin = start frequency as 10**Dmin
     Dmax = stop frequency as 10**Dmax
     mode = display mode: 'dB' magnitude, 'phase' in radians, or
            'groupdelay', all versus log frequency in Hz
     Npts = number of points to plot; defult is 1024
    fsize = figure size; defult is (6,4) inches

    Mark Wickert, January 2015
    """
    f = np.logspace(Dmin,Dmax,Npts)
    w,H = signal.freqs(b,a,2*np.pi*f)
    plt.figure(figsize=fsize)
    if mode.lower() == 'db':
        plt.semilogx(f,20*np.log10(np.abs(H)))
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Gain (dB)')
        plt.title('Frequency Response - Magnitude')

    elif mode.lower() == 'phase':
        plt.semilogx(f,np.angle(H))
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Phase (rad)')
        plt.title('Frequency Response - Phase')

    elif mode.lower() == 'groupdelay':
        """
        Notes
        -----

        See freqz_resp() for calculation details.
        """
        theta = np.unwrap(np.angle(H))
        # Since theta for an FIR filter is likely to have many pi phase
        # jumps too, we unwrap a second time 2*theta and divide by 2
        theta2 = np.unwrap(2*theta)/2.
        theta_dif = np.diff(theta2)
        f_diff = np.diff(f)
        Tg = -np.diff(theta2)/np.diff(w)
        max_Tg = np.max(Tg)
        #print(max_Tg)
        plt.semilogx(f[:-1],Tg)
        plt.ylim([0,1.2*max_Tg])
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Group Delay (s)')
        plt.title('Frequency Response - Group Delay')
    else:
        print('Error, mode must be "dB" or "phase or "groupdelay"')

Example: Discrete-Time Chebyshev Type I Bandpass Filter

import sk_dsp_comm.iir_design_helper as iird
import sk_dsp_comm.fir_design_helper as fird

b1,a1,sos1 = iird.IIR_bpf(200,250,300,350,0.1,60.0,1000,'butter')
b2,a2,sos2 = iird.IIR_bpf(200,250,300,350,0.1,60.0,1000,'cheby1')

IIR butter order = 16.
IIR cheby1 order = 12.

figure()
iird.freqz_resp_cas_list([sos1,sos2],'dB',1000)
ylim([-70,0])
grid();
figure()
iird.freqz_resp_cas_list([sos1,sos2],'groupdelay_t',1000)
grid();
figure()
iird.sos_zplane(sos2)

/home/docs/checkouts/readthedocs.org/user_builds/scikit-dsp-comm/envs/develop/lib/python3.5/site-packages/scikit_dsp_comm-1.1.0-py3.5.egg/sk_dsp_comm/iir_design_helper.py:346: RuntimeWarning: divide by zero encountered in log10
  plt.plot(f*fs,20*np.log10(np.abs(H)))
/home/docs/checkouts/readthedocs.org/user_builds/scikit-dsp-comm/envs/develop/lib/python3.5/site-packages/scikit_dsp_comm-1.1.0-py3.5.egg/sk_dsp_comm/iir_design_helper.py:379: RuntimeWarning: divide by zero encountered in log10
  idx = np.nonzero(np.ravel(20*np.log10(H[:-1]) < -400))[0]
/home/docs/checkouts/readthedocs.org/user_builds/scikit-dsp-comm/envs/develop/lib/python3.5/site-packages/scikit_dsp_comm-1.1.0-py3.5.egg/sk_dsp_comm/iir_design_helper.py:379: RuntimeWarning: invalid value encountered in multiply
  idx = np.nonzero(np.ravel(20*np.log10(H[:-1]) < -400))[0]

(12, 12)

<Figure size 432x288 with 0 Axes>

<Figure size 432x288 with 0 Axes>

<Figure size 432x288 with 0 Axes>

b,a = signal.cheby1(5,.1,2*array([250,300])/1000,btype='bandpass')

freqz_resp(b,a,mode='dB',fs=1000,fsize=(6,2))
grid()
ylim([-80,5]);
xlim([100,400]);
freqz_resp(b,a,mode='groupdelay_s',fs=1000,fsize=(6,2))
grid()
xlim([100,400]);

Example: Continuous-Time Bessel Bandpass Filter

bc,ac = signal.bessel(7,2*pi*array([10.0,50.0])*1e6,btype='bandpass',analog=True)

freqs_resp(bc,ac,6,9,mode='dB',fsize=(6,2))
grid()
ylim([-80,5]);
freqs_resp(bc,ac,6,9,mode='groupdelay',fsize=(6,2))
grid()

Second-Order Butterworth Lowpass Response

b3,a3 = signal.butter(3,2*pi*1,analog=True)
freqs_resp(b3,a3,-1,2,mode='dB',fsize=(6,2))
grid()
ylim([-80,5]);
freqs_resp(b3,a3,-1,2,mode='groupdelay',fsize=(6,2))
grid()

Obtaining the Step Response via Simulation

Time domain simulation of continuous time system can be performed using the signal.lsim() function. You have to make sure the time step is sufficiently small relative to the filter bandwidth.

t = arange(0,2,.0001)
xs = ss.step(t)
tout,ys,x_state = signal.lsim((b3,a3),xs,t)
plot(t,ys)
title(r'Third-Order Butterworth Step Response for $f_3 = 1$ Hz')
ylabel(r'Ste Response')
xlabel(r'Time (s)')
grid();